% Note: 
%   This function simulate a panel with entry and exit for the steady state 
%   distribution.
% Structure of the panel is a 3-D matrixL:Obs ID,Variable ID,Period T).
% In each period t, the cross-sectional data is structured as:
%
%   StateMat:       [b,a,zID,Age]
%   PolicyMat:      [bp,ap,c,l]
function SimuSample=SteadyState_Distribution_Simulation(PP,SS,N,T)
% N           =   5000;
% T           =   100;
%% Preliminary
% Preset the RNG 
rng(12345);

%--------------------------------------------------------------------------
% Entrant Distribution Parameters 
%--------------------------------------------------------------------------
% Idio. TFP (quadrature)
QW_z0       =   PP.ExoState.IdioInc.InitialDist;
% b
b0          =   0;
% a
a0          =   0;

%--------------------------------------------------------------------------
% Important Constants
%--------------------------------------------------------------------------
PolicyList  =   {'bp','ap','c','l','zl'};
N_Policy    =   length(PolicyList);
PolCoefMat  =   struct();
for ii=1:N_Policy
    vv              =   PolicyList{ii};
    TempPolicyMat   =   reshape(SS.Policy.V.(vv),SS.PolApp.V.MatSize_State);
    PolCoefMat.(vv) =   SS.PolApp.V.BasMat_0\TempPolicyMat;
end
%% Preset the History
Hist_State  =   zeros(N,4,T);
Hist_Policy =   zeros(N,4,T);
%--------------------------------------------------------------------------
% Exogenous Part of the History
%--------------------------------------------------------------------------
% zID and Age
Hist_zID_Age=   MarkovChainSimulation(full(PP.ExoState.IdioInc.TrMat),QW_z0,PP.XI, ...
                                      rand(N,T),rand(N,T));
% Entrant Distribution (b,a)
Hist_EntDist=   zeros(N,2,T);

%% Simulation of History
% Initialization
Temp_b_a    =   squeeze(Hist_EntDist(:,:,1));

% Simulation
for tt=1:T
    zID_Age     =   squeeze(Hist_zID_Age(:,:,tt));
    EntDist     =   squeeze(Hist_EntDist(:,:,tt));
    [PolicyMat,StateMat] ...
                =   SteadyState_Distribution_Simulation_OneStep( ...
                            PP,SS,zID_Age,EntDist,Temp_b_a,PolCoefMat);
    Temp_b_a    =   [PolicyMat(:,1),PolicyMat(:,2)];
    
    Hist_State(:,:,tt) ...
                =   StateMat;
    Hist_Policy(:,:,tt) ...
                =   PolicyMat;
end

%% Summarize the Ergodic Distribution
QW          =   struct();
QW_List     =   {'b','a'};
[QW.b.w,QW.b.Breaks] ...
            =   histcounts(StateMat(:,1),50,'Normalization','probability');
[QW.a.w,QW.a.Breaks] ...
            =   histcounts(StateMat(:,2),50,'Normalization','probability');
for ii=1:length(QW_List)
    vv          =   QW_List{ii};
    QW.(vv).w   =   [0;QW.(vv).w';0];
    QW.(vv).Breaks ...
                =   QW.(vv).Breaks';
    QW.(vv).Node=   [QW.(vv).Breaks(1); ...
                     (QW.(vv).Breaks(1:end-1)+QW.(vv).Breaks(2:end))/2;...
                     QW.(vv).Breaks(end)];
end

%% Combine the Results
SimuSample  =   struct('StateMat',StateMat,'PolicyMat',PolicyMat, ...
                       'Hist_State',Hist_State,'Hist_Policy',Hist_Policy,...
                       'QW',QW);
                   
end

%   StateMat:       [b,a,zID,Age]
%   PolicyMat:      [bp,ap,c,l]
function [PolicyMat,StateMat] ...
        =SteadyState_Distribution_Simulation_OneStep(PP,SS,zID_Age,EntDist,Temp_b_a,PolCoefMat)
%% Preliminary
N_State     =   size(zID_Age,1);
%% One-Step Transition
% Initialization
zID         =   zID_Age(:,1);
b           =   Temp_b_a(:,1);
a           =   Temp_b_a(:,2);
% Exogenous Exit
TempInd     =   ( zID_Age(:,2)==0 );
b(TempInd)  =   EntDist(TempInd,1);
a(TempInd)  =   EntDist(TempInd,2);

% Portfolio Choice and Consumption
PolicyState =   [b,a,zID];
bp          =   zeros(N_State,1);
ap          =   zeros(N_State,1);
c           =   zeros(N_State,1);
l           =   zeros(N_State,1);
for ii_z=1:PP.ExoState.IdioInc.N
    TempInd     =   (PolicyState(:,3)==ii_z);
    if any(TempInd)
        TempBasMat      =   funbas(SS.PolApp.V.Space,PolicyState(TempInd,1:2),[0,0]);
        bp(TempInd)     =   TempBasMat*PolCoefMat.bp(:,ii_z);
        ap(TempInd)     =   TempBasMat*PolCoefMat.ap(:,ii_z);
        c(TempInd)      =   TempBasMat*PolCoefMat.c(:,ii_z);
        l(TempInd)      =   TempBasMat*PolCoefMat.l(:,ii_z);
    end
end

%% Collect the PolicyMat and LagStateMatp
PolicyMat   =   [bp,ap,c,l];
StateMat    =   [b,a,zID_Age];

end